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Summary 


T. J. R. Hughes and W. K. Liu described a new type of implicit explicit 
algorithm for the solution of transient problems in structural dynamics. The method 
involved dividing the finite elements into implicit and explicit groups while 
automatically satisfying the interface conditions. It is the purpose of this re- 
port to apply this algorithm to the solution of the linear, transient, two-dimen- 
sional heat equation subject to an initial condition derived from the solution of a 
steady state problem over an L-shaped region made up of a good conductor and an 
insulating material. 

Using the IIT/PRIME computer with virtual memory, a FORTRAN computer program 
code was developed to make accuracy, stability, and cost comparisons among the fully 
explicit Euler, the Hughes-Liu, and the fully implicit Crank-Nichol son algorithms. 
This report illustrates the Hughes-Liu claim that the explicit group governs the 
stability of the entire region while maintaining the unconditional stability of the 
implicit group. 


Introduction 

Transient heat flow in solids can be described mathematically in terms of a 
parabolic differential equation. For certain very simple boundary conditions, this 
equation may be solved analytically, but when complex problems are considered, this 
is usually not possible. In that case the transient heat transfer problem may be 
descretized by finite elements and a weak solution is obtained at the nodes of the 
mesh. The weak solution is adequate provided it converges. 

The truncation error, (the difference between the exact solution of the 
transient problem, and the exact solution of the discretized problem) is due to the 
finite distance between mesh points, the approximation of the time derivatives by 
finite differences, and the choice of the integration scheme. To find conditions 
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under which the truncation error goes to zero, is the problem of convergence [1]. 

By a well-known result, stability and consistency imply convergence, so that sta- 
bility can be a computationally crucial question in such problems. 

There are two major types of algorithms for the solution of transient structural 
heat transfer problems: implicit and explicit. In implicit algorithms the nodal 

temperatures at time t + At are expressed as functions of temperatures as time t and 
t + At, while in explicit algorithms they are expressed solely as functions of 
temperatures at time t. Explicit algorithms require much less computation per time 
step, but the time step size is limited (often severely) by the stability considera- 
tion. These stability limitations are most severe for good conducting materials. [2]. 

The present work considers the solutions of the transient heat transfer problem 
in an L-shaped region, one part of which consists of a very good conductor, and the 
other part is made up of insulation. In this type of situation researchers have 
begun to consider the possibility of solving the problem by a mixed implicit-explicit 
algorithm. Belytschko, and Mullen [3] have an implicit method that works on the 
following idea of strong coupling: The explicit group is integrated first, and the 

information then moves between adjacent nodes during a particular time step. This 
then is used as a boundary condition to integrate the implicit group of elements, 
and information in the implicit group travels across the entire mesh. In order to 
accomplish this, it is necessary to partition the element into implicit, explicit, 
and interface groups, while partitioning the nodes into implicit and explicit 
groups [3]. 

Hughes and Liu [4,5] have a simpler implicit-explicit algorithm where the 
elements are subdivided into implicit and explicit groups, and no provision needs to 
be made for interface elements since these are automatically taken care of during 
the assembly process. It is the purpose of this report to apply the Hughes-Liu 
algorithm to the heat equation, and to make appropriate analysis of both stability 
and accuracy. The performance of the Hughes-Liu algorithm is then compared with 



the implicit Crank-Nicol son and the explicit Euler methods. These three methods 
are also compared from the point of view of cost in terms of the CPU time {all cal- 
culations were performed on the I IT prime 400 computer). 

Problem Description 


Heat flow in solids in the absence of internal heat generation is described by 
the equation: 


V • (kvu) = cu^ 
k = thermal conductivity 
c = specific heat capacity 



( K) 


0 ) 


y = temperature 
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t = subscript on u denoting differentiation 
with respect to time in seconds 

Equation (1) is solved herein by the finite element method over the L-shaped 
region composed of two materials. The geometry and boundary conditions are indicated 


in Fig. 1 . 
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Figure 1 

The L-shaped Region 
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The spatial discretization of Figure (1) is accomplished by a "crude" mesh of 
twelve elements (twenty one nodes) of dimensions (1/16) by (1/16), and (1/16) by 
(7/16), and by a "fine" mesh of two hundred and forty elements (three hundred and 
five nodes) of dimension (1/32) by (1/32). At the element level, both k and c are 
constant, so that at each point of the region, equation (1) may be simplified when 


expressed in terms of the diffusion coefficient a given by 

.2 


k m‘ 

oi = zr 


c sec 


( 2 ) 


so that 


aV u = u. 


(3) 


In the present work the initial condition: 


u(0) = U (4) 

is obtained by solving the steady state problem: 

av^u = 0 (5) 

with the same boundary conditions. By means of the Galerkin Weighted Residual 
method, for example, equation (3) is reduced to the equivalent matrix equation: 

Cu + Ku = 0 (6) 


where : 


C = global heat capacity matrix 
K = global thermal conductivity matrix 
u = vector of unknown nodal temperatures 

and the dot denotes time differentiation. Also the steady state problem (5) is re- 
duced to the matrix equation: 


Ku = 0 


(7) 


The solution of Eq. 7 subject to the boundary condition in Fig. 1 gives the initial 
condition used in Eq. 4. 
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Explicit and Implicit Algorithms 
The following notation is used in the solution algorithm: 


d = u 


(8) 

(9) 

00 ) 


In equations (9) and (10) d^ and u^ are the numerical approximations to the exact 


value at time 


where 


= nAt (11 ) 


t^ = n^*^ timestep 
At = time increment 

The explicit algorithm used is the Euler method which consists of the first two 
terms of the Taylor series. Thus it is a first order method and is given by: 


"n+l = “n ^ 


(12) 


The stability of this method may be characterized in terms of the Fourier number 
which is defined as the ratio of the rate of heat conduction to the rate of heat 
storage. In physical terms it gives an indication as to how fast the temperature of 
the material changes due to a heat input [7]. If Ax, and Ay are the minimum element 
lengths in the model of the region along the x, and y directions, respectively, then 


the Fourier number, F^, is defined by: 


where 


F = aAt 
0 


a = 


(ax)^ + (Ay)^ 


(13) 


(14) 


In the present problem we are dealing with an inhomogeneous material region 
consisting of a good conductor and an insulator. We have found that the size of At 
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is limited by the higher diffusion ratio of the good conductor. Furthermore, At 
must be chosen small enough to meet the Courant stability criterion: 

< % (15) 

The implicit algorithm used is the Crank-Nichol son method. It is a second order 
method given by: 

which is unconditionally stable. In the numerical experiments the solution obtained 
with the implicit method with the time step of the explicit method is considered to 
be standard (baseline) for accuracy considerations. 


The Hughes-Liu Algorithm 


The elements of the mesh are divided into implicit and explicit groups. 
Hughes-Liu algorithm replaces equation (6) by: 


The 


Cd + K^u + K^u = 0 
n+1 n+1 n+1 


(17) 


where 


I indicates implicit group 
E indicates explicit group 


'\j 


The predictor value u for u is 


u„ + %Atd^ 

n+1 n n 


(18) 


which is just the central difference formula for d . The corrector equation is 


'n+1 " "n+1 * '^**^'^n+l 


(19) 


Reduction of Hughes-Liu Algorithm to a Single Equation 

The predictor equation (18) is solved for d^ and the corrector equation (19) is 
evaluated at the n^'^ step. Then with d^ eliminated from equation (19) the following 
predictor-corrector relation is obtained: 
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u = %(u + u , 1 ) 

n n n+r 


(20) 


With equation (20) the equivalent split operator problem (17), (18), and (19) is 
reduced to the single equation: 


Vl = (21) 

where 

G(At) = (I - %AtA^)"^(l + %AtA^ + AtA^) (22) 

wi th 

A^ = - (23) 

E -IF 

A = - C K (24) 

I E 

For the fully explicit case (K = 0, K = K) equation (25) simplifies to: 

G(At) = I + AtA (26) 

where 

A = A^ (27) 

Thus the fully explicit equation (26) is the first order Euler method. For the 
fully implicit case (K^ = K, = 0) equation (22) simplifies to: 

G(At) = (I - %AtA)'^(I + %AtA) (28) 

where 

A = A^ (29) 

and this is the second order Crank-Nichol son method. 


Stability Analysis for Hughes-Liu Algorithm 

Equation (21) may be put into the form: 

u^^^ = (I - AtB)u^ (30) 

(See Appendix for a one-dimensional example). Equation (30) represents a contraction 
map if for the largest eigenvalue of B we have: 

< 1 


1 - AAt 


(31) 
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Using an energy type norm of the discrete solution, Hughes and Liu show for a 
dynamics problem that the largest eigenvalue of B is the same as if only the explicit 
region was analyzed [4]. For heat conduction this is demonstrated in the Appendix 
for the one-dimensional case. Because the explicit region is chosen to be the in- 
sulator, the stability circle for the Hughes-Liu method is considerably larger than 
the unit circle centered at -1 which is the stability circle for the explicit Euler 
method. This is shown in Figure (2). The fully implicit Crank-Nichol son method is 
unconditionally stable, so it has a stability region consisting of the entire left 
half of the complex plane. (See Figure (3). 



Figure (2) 

Stability regions for the Euler and Hughes-Liu methods 
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Im(XAt) 


Re(AAt) 


Figure (3) 

Stability region for the Crank-Nichol son method 


Steady State Solution 

Figure (4) shows the coarse mesh discretization of the L-shaped region and the 
steady state distribution of the initial temperature which is input to the transient 
problem. Material properties are given in Table 1 and the value of the Fourier con- 
stant (Eq. 14) in Table 2. The temperature distribution is tabulated for diffusion 
ratios of one, fifty, and one hundred. While the normal derivative is zero along 
the edge between nodes 14 and 15, it has the approximate value of -1.1 (K/m) along 
the edge between nodes 15 and 10, resulting in a sharp discontinuity at node 15. 
Careful examination of Figure (8.1) shows the while the temperature along the 
essential boundary nodes (5-10-15) decreases lOOK — >■ 70.71 K — OK, the temperature 
along the opposite end of elements 4 and 8 along nodes (4-9-14) increases 
36.52K — ^ 43.72K — 62.23K. This phenomenon of "reflection" is the direct result 
of the normal derivative discontinuity at the corner node 15 and it has been re- 
solved by mesh refinement. In particular, with a finer mesh, where elements 4 and 
8 of the coarse mesh have been replaced by the 56 elements, at the position cor- 
responding to nodes (4-9-14) the temperature became constant at 47.70K. 




: DISTRIBUTION FOR COARSE MESH 


Figure 4 


Di ffusion 
Ratio 

Nodal Numbers* 

(16) . (17) • (18) 

1 

13.74 

13.76 

13.98 

50 

28.26 

28.31 

29.54 

100 

28.59 

28.64 

29.91 


*Table continued below 


(14) 

[7] 

(9) 

[8] 

[3] 

(4) 

[4] 


(15)0K 


(12) 

(3) 

(8) 

(13) 

(4) 

(9) 

27.81 

58.40 

59.12 

36.30 

62.76 

63.34 



33.87 

60.39 

60.96 

26.92 

53.29 

53.84 

36.52 

49.76 

50.04 

43.72 

56.94 

57.22 
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Experiment 1: Accuracy of Fully Explicit vs. Hughes-Liu Algorithm 

The purpose of this experiment was to compare the accuracy of the Hughes-Liu 
method and fully explicit method for the same step size. The fully implicit method 
which is of order two (compared to order one for the fully explicit) is taken as the 
standard. 

For a fixed time step the conductor has a Fourier number two orders of magnitude 
larger than that of the insulation. For this reason implicit integration is re- 
quired over the conductor while explicit integration over the insulation seems 
sufficient. The time step must be chosen small enough to meet Courant's stability 

criterion (F^ < for the conducting material. In particular, in case of the fine 

_2 

mesh, the conducting material has a = 3.2 x 10 and for At = 7.81 seconds the value 
of Fq is \ as shown in Table 1. 

As shown in Table 3, after 77 steps, (elapsed time is 602 seconds), the maximum 
error of 0.21 K occurred at nodes 26, 92, and 258 of the fine mesh for the fully ex- 
plicit method relative to the standard. This is a negligible error occurring over 
the conductor near its essential boundary which illustrates that the Courant con- 
dition is sufficient for stability of explicit time integration. The Hughes-Liu 
method showed the same accuracy as the standard but had the advantage in CPU time: 

983 sec. for the Hughes-Liu method vs. 1327 sec. for the fully implicit. 

While the fully implicit method is unconditionally stable convergence is not 
guaranteed. In a run with At = 600 sec., negative temperature values appeared near 
the essential boundary of the conductor when calculated for the first and only step. 
The reason for this is the presence of large truncation error (the Fourier number 
is 19.2). 


Experiment 2: Determination of Maximum Time Step 


The purpose of this experiment was to determine the maximum time step for each 



method to achieve accuracy within one degree of the temperature. The baseline for 
this coincided with the standard at the time step of 7.8 sec., and the final time 
was 78 sec. The material properties for the insulator were varied to achieve 

4 

various diffusion ratios by fixing heat capacity at 9.4 x 10 and adjusting thermal 
conductivity as shown in Table 4. Using the standard. Table 5 shows that varying 
the diffusion ratio has very little effect on the baseline temperature, or on the 
variation from the baseline by the three methods. 

The data presented in Table 5 is at nodes 31, 64, 97, 130, and 163, located 
1/16 meters from the essential boundary of the conductor. These nodes were selected 
because the maximum temperature variation within the one degree limit had been ob- 
served there. It can be seen in Table 5 that the fully explicit method is within 
one degree of baseline using a time step of 7.1 sec. This corresponds to = 0.2273 
just 9.09% below the Courant value. In case of the fully implicit method the cor- 
responding time step is 26.04* sec., which is 3.33 times larger than the fully ex- 
plicit time step. In the next two experiments with larger material properties this 
ratio was increased to 9.00. The Hughes-Liu method was compared to the baseline only 
at the diffusion ratio of fifty. It was found to be identical to the fylly implicit 
results shown in Table 5. The CPU times for the three methods are shown in Table 6. 
It is seen that the Hughes-Liu method is the least expensive of the three methods. 

Experiment 3: Effect of the Diffusion Ratio 

In this experiment, the thermal conductivity for the good conductor was in- 
creased from its value in Experiment 2 to 1500 ( Watt/m- K)), the time step was 
correspondingly reduced to 0.78 sec., and the thermal conductivity was increased by 
a factor of ten over those listed in Table 4. The data is presented at the nodes 
32, 65, 131, 98, and 164, located 1/32 meter from the essential boundary of the 

* Note that the choice of a time step is constrained by the requirement that the 
final time is an integer multiple of the time step. 
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conductor, slightly closer than in Experiment 2, because there is where the maximum 
variation has shifted as the result of the above changes while remaining below one- 
degree of baseline. Table 7 shows that the effect of changing the diffusion ratio 
is negligible. It is seen that baseline temperatures are identical at each node. In 
the fully implicit case of Table 7, a difference of 0.01 K shows up at nodes 65 and 
131 as compared to Table 5 with nodes 31, 64, and 130. 

The fully explicit method is within one degree of baseline at t = 0.79 sec., 
corresponding to = 0.2525, which is just 1.01% above the Courant value of 0.25. 
The fully implicit method remains within one degree of baseline at the same nodes 
when the time step equals 7 sec. The ratio between the fully implicit and fully 
explicit time steps is 9.0. 

The Hughes-Liu method was compared to the baseline only at the diffusion ratio 
of fifty, and the temperature values were found to be identical to the values ob- 
tained in the fully implicit case. The CPU times for the three methods are shown 
in Table 8 and show that while the CPU time for the fully implicit and Hughes-Liu 
method doubled, the CPU time for the fully explicit method increased by a factor 
of 7.65, making it much more expensive than either of the other two methods. Since 
the Hughes-Liu method is as accurate as the fully implicit method, the Hughes-Liu 
is the superior method for the present problem. The increased efficiency of the 
implicit method and the Hughes-Liu method over the explicit method is due to the 
increased diffusivity . Because of the higher value of the diffusivity the tem- 
peratures decay faster from the initial state to the final state. As a result, a 
larger share of the temperature history is in a region where there is little 
temperature variation as compared to the previous experiment. 

Experiment 4: Effect on Implicit Algorithm of Variable Time Steps 

The purpose of this experiment was to find the effect of a variable time step 
on the implicit method. Because there is applied heat loading, accuracy require- 
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ments dictate small time steps at the beginning of the temp history. Later, as the 
temperatures become more uniform, the time step may be increased. For the implicit 
and Hughes-Liu methods the penalty for variable time step is the need to refac- 
torize the Jacobian. We thus used the crude mesh in order to minimize the cost of 
refactorization. Table 9 summarizes the results obtained for this experiment. It 
shows the effect on accuracy of using six integrating steps of variable size. Node 
14 was selected because that is where the maximum temperature variation from the 
baseline was observed. It is seen from Table 9 that by taking small time steps 
initially, the deviation of -10.66 K at fifty sec. is drastically reduced to -5.15 K 
at 53 sec. and temperature values more accurate. However, the three refactorizations 
required to achieve this increased accuracy costs more in terms of the CPU time. 

More specifically, while a single factorization run takes 3 sec., two additional 
refactorizations raise the CPU time to 5 sec. 

The number of steps to reach 300 sec. was increased from six to eight. The 
results are given in Table 10. Table 10 shows that by taking several small time 
steps of three seconds initially, the huge variation of 9.30°K at 37.5 sec. is 
reduced to -4.42 K at 40.5 sec., -2.35 K at 43.5 sec., and -1.32 K at 46.5 sec. 

While the one factorization run costs 3 sec. the (2-4-2) run costs 5 sec., and the 
(1-6-1) and (3-2-3) runs cost 6 sec. The effect of taking eight instead of six 
steps was that the initial oscillations are slightly smaller. This is simply due 
to the fact that in this case the maximum step size is 37.5 sec. as opposed to the 
50 sec, step size in Table 9. Tables 9 & 10 both show that for parabolic problems 
the initial oscillations do not affect the final temperature because they die 
out sooner. 

Concl usions 

This report describes the application of the Hughes-Liu algorithm to the solu- 
tion of linear, transient heat equations, subject to an initial condition. The 
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technique was applied to an L-shaped region made up of a good conductor and an 
insulating material. The report has illustrated that the explicit time step 
determined by the largest eigenvalue of the explicit region governs the stability 
of the entire region while unconditional stability is achieved for the implicit 
region. Numerical experiments have verified that this stable behavior leads to 
substantial computational savings in the two dimensional problem. 

Numerical examples have illustrated that when the diffusivity of the good 
conductor is high enough to warrant the use of a fully implicit integrator, while 
the poor conductor can be stably treated by explicit integration, the Hughes-Liu 
algorithm can lead to substantial computational savings. Caution should be ob- 
served in interpreting these findings because of the simple, rapidly decaying 
character of the test problem. Nevertheless, the implicit/explicit integration 

methods show promise as tools in the numerical solution of transient thermal 
analysis problems. 
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APPENDIX 

One Dimensional Stability Analysis for the 
Time Integration Methods 

If x-j , ^2 are sampling points; w-j , Wg are weights for two-point Gaussian 
quadrature rule over the closed interval from zero to one, then 
1 

F(s)ds = w-|F(x^) + W 2 F(x 2 ) (A1 ) 


I 

/ 


0 


If, further, we require that we have equal weights w, and symmetrically dis- 
tributed sampling points x, and 1 - x, then 
1 

F(s)ds = w[F(x) + F(1 - x)] (A2) 

0 

Being a two-point formula, the quadrature must be exact for polynomials of 
degree three, which leads to 
1 


I 

/ 


/ 


0 


= ^fFfi - ^ ) + Ffl + ~ 


F(s)ds = %[F(^ - -| ) + F( 


2 6 


)] 


(A3) 


I 

/ 


Finally, by moving the sampling points to the endpoints of the interval, 
we obtain the "lumped" quadrature rule [6] 

1 

F{s)ds %[F(0) + F{1)] (M) 

0 

The special notation used is listed below: 

k = Thermal conductivity for conductor 

a-|= Insulator to conductor conductivity ratio (0 < a-j < 1) 
c = Heat capacity for conductor 
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0 ^= Insulator to conductor capacity ratio (0 < Q 2 < 1 ) 
h = Element length 
X = Spatial coordinate 
t = Elapsed time 
u = Temperature 

I,E= Superscripts used to identify Implicitland Explicit groups 

Special Notation 


The mathematical 
k(x) = 


formulation : 
k 0 < X < h 

a^k h < X < 2h 


c(x) = 
k(x)u„^ = c(x)u 


I02C 


0 < X < h 
h < X < 2h 


XX 


0 < X < 2h, t > 0 


(A6) 


(A7) 

(A8) 


Subject to: 

Initial Condition 

u(x,0) = U(x) 0 < X < 2h (A9) 

BIC Conditions: 

(i) Essential: u(2h,t) =0 t > 0 (AlO) 

(ii) Natural: u (0,t) =0 t > 0 (All) 

A 

(iii) Compatibility: u^(h,t) = u^{h,t) t > 0 (A12) 

T E 

(iv) Interface: u^(h,t) = a,u (h,t) t > 0 (A13) 

A I A 

Solution Space: Hilbertian Sobolev Manifold 


(i) Trial function: hJ = |u e L^[0,2h]:u, u^ e BIc} (A14) 

(ii) Test function: hJ* = {v e L^[0,2h] :v,v^ e BIc} (A15) 

= hJ for Galerkin formulation 



I 


19 


Korn Inequality: 
where a , 3 are constants 


a u ^ < a(u,u)®< 3 u -| (A16) 


2 h 


(u^ + U^)dx 


(A17) 


and a(u,v) e / k(x)u vdx 

XX 


h 

/ 




induces the energy norm equivalent to 


If 


(A18) 


f(v) E 


c(x)u^vdx 


(A19) 


then the Galerkin method consists of solving for the residual 

R E a(u,v) - f(v) = 0 (A20) 

Using isoparametric transformation between global and local coordinates given 


by 

X = x(s) = hs + x^ 0 < s ^ 1 (A21) 

and assuming linear variation of temperature function u» u and x may be written in 


terms of shape functions 

(f)-| (s) =1 - s 
<l> 2 (s) = s 

so that 

a(u ,v ) = - V K u 

0,/ e^ T^e 
f(v ) = V C u 


(A22) 


(A23) 

(A24) 


20 


where 


.e _ r 
' T 


-1 


C® = he® 


-1 

1 


(A25) 


f (l--s)^ds f s(l-s)ds 
0 

/ 


0 

s(l-s)ds / s^ds 
0 


(A26) 


Applying integration formula (A4) equation (A26) becomes 

e 


^e _ he 


-1 

0 


(A27) 


Then separating out the formulas for the implicit and explicit groups, we obtain 


.1 _ k 
" h 


1 

-1 


o^k 

^ -IT 


pi _ he 
L 2 


^ o 


1 

-1 

1 

0 


- 1 " 

1_ 

-1 

1 

0" 

1 


which leads to the assembled matrices: 


K - - 
h 


r il£ 
^ “ 2 


1 

-1 

0 

1 

0 

0 


-1 

1+a 

-^1 

0 

1 +a, 

i 

0 


0 

1 '^1 


'1 


(A28) 

(A29) 

(A30) 

(A31) 


(A32) 


(A33) 
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The essential boundary condition = 0 forces the elimination of the 3rd row and 
3rd column from K, and C resulting in: 


K 


b 



-1 

1 +a-| 


C 


b 



0 

1 +02 


(A34) 


(A35) 


and K|^ can be split into implicit and explicit parts: 



(A36) 



whose eigenvalues are given by: 



Employing simultaneous diagonal ization on 
into the form 

Un+, = (I - BAt)u„ 


(A37) 

(A38) 

(A39) 

(A40) 

(A41) 

the Hughes-Liu equation, it may be put 


(A42) 
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and this same form may be achieved in both partially and fully explicit cases where 


- /TT 


2 nE 


1 + a« + q At 


/I + 

2q^ + q^ 

1 + + q^At 


B = 


- A 


In Hughes-Liu Equation 
0 

-Ao 


(A43) 


In partially Explicit Case (K^ = 0) 


2k 

h^c 


-1 


1 + a. 


-1 


1 + a-| 
1 + a. 


In Fully Explicit Case (kJ = 0, = Kj^) 

With corresponding eigenvalues: 

r - /p2 - s In Hughes-Liu Case 

E 


^ -^1 


In Partially Explicit Case 


2[r* - /(^*j2 _ 5 * In Fully Explicit Case 


^2 1 "^2 


r + /r2 - s 
E 


In Hughes-Liu Case 
In Partially Explicit Case 


2[r* + /(^*)2 _ 2 * In Fully Explicit Case 

h c 


(A45) 
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Where : 


I 

q 

= (2 + 0 ^) 2 
h c 


(A46) 

E 

k 


(A47) 

q 

^ ^ °2 h^c 


V*i — 

q^q^At + (2 + 

02)q^ + 2q^ 

(A48) 

1 “ 

2(q^At + 

1 + 1 ^ 2 ) 


2q q 


(A49) 

S " 

q^At +1+02 


r* 

1+0, 

) 

(A50) 

s* 

1 + O 2 


(A51) 

The stability condition corresponding to equation (A42) is given by 

|i 

- AtA| < 1 


(A52) 

and since in all cases 




^1 

< X 2 


(A53) 

is equivalent to 




0 5 

At < ^ 

^2 


(A54) 

so that 

' 1 ^ “2 h^c 

o-| k 

in Hughes-Liu Case 


0 < At < ^ 

1 1 " ^2 h^c 

/ 0 ^ k 

In partially Explicit 

Case 


i 2 + + o^ - 

■'(°,-a2)^+4(l+a2) h^c 



' 

2o^ k 



\ In Fully Explicit Case 

(A55) 
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This illustrates the energy theorem of Hughes and Liu [4] which states that when- 
ever the partially explicit case is stable over the explicit region, the Hughes-Liu 
method will be stable over the entire inhomogeneous region. 

Finally, in the fully implicit case (K^ = 0, kJ = Kj^) we obtain the equation 

(I = (I - ^B)Ur, (A56) 

where B is the matrix obtained in the fully explicit case. 

This shows that the stability criterion is 

5 1 (A57) 

so that At > 0 unconditionally. 
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Table 1 


Material Properties of L-shape Configuration 



Thermal 

Heat 

Di ffusion 

Fourier 

Material 

Conductivity 

Capacity 

Coefficient 

Number, F^, Eq. 13 

Conductor 

150 

2,4 X 10^ 

6.25 X 10'^ 

0.25 

Insul ation 

0.06 

9.4 X 10^ 

6.38 X 10'^ 

0.00255 


Table 2 

Value of Fourier Constant a (Eq. 14) 











Table 3 

Comparison of Accuracy of the 

Fully Implicit, Fully Explicit, and Hughes-Liu Methods at Selected Nodes 


Time Step =7.8 sec. 
Number of Steps = 77 
Duration - 602 sec. 


Node 

Number 

Implicit 

Expl icit 

Error 

Hughes 

Liu 

Error 

1 

63 

63 


63 

0.00 

2 

63 

63 

0.00 

63 

0.00 

3 

63 

63 

0.01 

63 

0.09 

26 

36.13 

36 

0.21 

36.13 

0.00 

92 

36.13 

36 

0.21 

36.13 

0.00 

258 

20.19 

20.19 

0.00 

20.19 

0.00 

CPU 

TIME 

1327 
sec . 

507 
sec . 


983 
sec . 

TOTAL 
2939 
sec . 





























28 


Table 4 

Insulator Thermal Conductivity for Various Diffusion Ratios: 


Di ffusion 

Thermal 

Rati 0 

Conductivity 

1 

5.875 

50 

0.1175 

100 

0.05875 


Table 5 

Temperature Variation from Baseline for 

Various Diffusion Ratios (At chosen to achieve 
errors below IK) 


■i 



Temperature 

Variation from 




Baseline 



Basel ine 



Ml 

Di ffusion 

Temperature(K) 

Fully Implicit* 

Fully Explicit 


Ratio 

At = 7.8 sec. 


At = 7.1 sec. 


1 

29.71 

0.78 

-0.74 

31 

50 

29.87 

0.78 

-0.73 


100 

29.87 

0.79 

-0.73 


1 


0.71 

-0.71 

64 

50 


0.70 

-0.72 


100 

29.84 

0.71 

-0.72 


1 

29.58 

0.48 

-0.67 

97 


29.75 

0.48 

-0.68 


100 

29.75 

0.48 

-0.67 


1 

29.50 

0.22 

-0.63 

130 


29.66 

0.23 

-0.63 


100 

29.67 

0.22 

-0.63 


1 

29.46 

0.09 

-0.61 

163 

50 

29.63 

0.09 

-0.62 


100 

29.63 

0.09 

-0.62 


*The Hughes-Liu method produced identical results for diffusion ratio of 50 
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Table 6 

CPU Times for various algorithms diffusion ratio = 50 

(Experiment 2) 


Integration 

CPU 

Method 

TIME 

Fully Implicit 

127 sec. 

Fully Explicit 

33 sec. 

Hughes-Liu 

77 sec. 


Table 7 

Temperature Variation from Baseline for 
Various Diffusion Ratios (Experiment 3 
increased conductivity, At chosen to 
achieve errors below IK) 


Node 

Di f fusion 
Rati 0 

Basel i ne 
Temperature (K) 
At = 7.81 sec. 

Temperature Variation from 
Basel ine 

Fully Implicit* 
At = 7.1 sec. 

Fully Explicit 
At = 0.79 sec. 


1 

4.92 

-0.59 

-0.10 

32 

50 

5.01 

-0.60 

-0.11 


100 

5.01 

-0.60 

-0.11 


1 

4.92 

-0.55 

-0.11 

65 

50 

5.01 

-0.56 

-0.12 


100 

5.01 

-0.55 

-0.12 


1 

4.92 

-0.43 

-0.11 

98 

50 

5.01 

-0.44 

-0.12 


100 

5.01 

-0.44 

-0.11 


1 

4.92 

-0.20 

-0.10 

13T 

50 

5.01 

-0.21 

-0.11 


100 

5.01 

-0.20 

-0.11 


1 

4.92 

-0.14 

-0.11 

163 

50 

5.01 

-0.15 

-0.12 


100 

5.01 

-0.15 

-0.12 


*The Hughes-Liu method produced identical results for diffusion ratio of 50 



30 


Table 8 

CPU Times for Experiment 3 at diffusion ratio = 50 


Integration 

CPU 

Method 

TIME 

Fully Implicit 

245 sec. 

Fully Explicit 

635 sec. 

Hughes-Liu 

158 sec. 


Table 9 

Effects of variable time step on error in 
temperature at Node 14 (six time steps, implicit method) 


Elapsed 

Time 

sec. 

Basel ine 
at 

t=l sec. 

6 steps at At=50 sec. 

1 steps at At=3 sec. 

4 steps at At=50 sec. 
1 steps at At=50 sec. 

2 steps at At=3 sec. 

2 steps at At=50 sec. 
2 steps at At=97 sec. 

3 

67.89 


-0.27 

-0.27 

6 

63.67 



-0.19 

50 

57.89 

-10.66 



53 

51 .45 


-5.15 


56 

51 .03 



-2.78 

100 

45.56 

7.25 



103 

45.23 


3.24 


106 

44.90 



1 .64 

150 

40.66 

- 5.15 



153 



-2.30 


200 


3.52 



203 



1 .42 

-1.56 

250 

33.47 

- 2.55 



300 

30.69 

1 .87 

-1.33 

0.97 














Table 10 


Effect of Variable time steps on temperature error (node 14) 


Elapsed 
T i me 
in 
Sec. 

At=.05 sec.' 

8 steps At=37.5 sec. 

1 step At=3 sec. 

6 steps At=37.5 sec. 
1 step At=72 sec. 

2 steps At=3 sec. 

4 steps At=37.5 sec. 
2 steps At=72 sec. 

3 steps At=3 sec. 

2 steps At=37.5 sec. 
: 3 steps At=72 sec. 

3.0 

67.91 


-0.29 

-0.29 

-0.29 

6.0 

63.68 



-0.20 

-0.20 

9.0 

61.20 




-0.14 

37.5 

53.79 

-9.30 




40.5 

53.31 


-4.42 



43.5 

52.85 



-2.35 


46.5 

52.40 ! 




-1 .32 

75.0 

48.50 ! 

5.42 




78.0 

48.13 1 


2.45 



81.0 

47.76 



1 .22 


84.0 

47.40 




0.64 

112.5 

44.22 

-3.60 




115.5 

43.91 


-1 .52 



118.5 

43.61 



-0.77 


150.0 

40.66 

2.26 




153.0 

40.40 


0.83 



156.0 

40.14 



0.36 

0.61 

187.5 

37.64 

-1.60 




190.5 

37.41 


-0.59 



225.0 

35.03 

1 .03 




228.0 

34.83 


0.32 

-0.38 

0.26 

300.0 

30.69 

0.52 

-0.33 

0.18 

-0.35 
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